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SUMMARY 

The singularity problem associated with a radially continuous Maxwell viscoelastic 
structure is investigated. A special tool called the isolation function is developed. 
Results calculated using the isolation function show that the discrete model assumption 
is no longer valid when the viscoelastic parameter becomes a continuous function of 
radius. Continuous variations in the upper mantle viscoelastic parameter are especially 
powerful in destroying the mode-like structures. The contribution to the load Love 
numbers of the singularities is sensitive to the convexity of the viscoelastic parameter 
models. The difference between the vertical response and the horizontal response found 
in layered viscoelastic parameter models remains with continuous models. 
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1 INTRODUCTION 

In his remarkable monograph Theory of Viscoelasticity, an 
Introduction’ Christensen (1971) warned that the extension of 
methods developed for homogeneous materials to inhomo- 
geneous matenais may be difficult, if not impossible. He did 
not, however, detail how and when the difficulties would occur. 
Three years later, in studying the viscoelastic response of the 
mantle to deglaciation, Peltier (1974) published perhaps the 
first such extension for a Maxwell viscoelastic rheology. In his 
paper, Peltier (1974) assumed that, for an inhomogeneous 
viscoelastic parameter t](r)/fi{r), which is a piecewise continuous 
function of radius, r, the relaxation of the Maxwell mantle is 
in the form of a set of discrete exponentially decaying time 
functions exp(s ; -r), s } < 0, j = 1, 2, 3, • ■ ■ (modes). This assump- 
tion is warranted in an exact sense when r}(r)jjx{r) is formed by 
homogeneous layers (e.g. Wu & Peltier 1982). Each viscosity 
jump introduces two ‘viscosity’ modes (see below) and each 
density jump causes one ‘buoyancy’ mode (Han & Wahr 1995). 
The secular equation for these modes can be written as 

21 V - Si 

I v* = o, 

n = 0 

where N and M are the numbers of viscosity jumps and density 
jumps, respectively. If we follow the spirit of using more and 
more layers to approximate a continuous structure, we encoun- 
ter a theoretical crisis at its limit, namely, the solution of the 
secular equation becomes meaningless when 2JV + M tends 
to infinity. Practically, it is very difficult to denve analytic 
expressions for the coefficients ct„ even for fairly small V, say 
<V = 10. The solution of the resulting high-order algebraic 
equation is also a big problem. Hence, most people turn to 


the Runge-Kutta propagation technique, commonly called the 
normal mode method (in this paper, the term ‘normal mode 
method’ is only used in relation to the propagation technique), 
to search for the discrete modes (e.g. Wu and Peltier 1982; 
Peltier 1985; Han & Wahr 1995). In this way, the ‘theoretical 
crisis’ mentioned above shifts to a singularity problem. 

To be specific, let 


H(r) 

, = max 

Hir) 

n(r)_ 

y max 

w). 


Then, for a fixed Laplace transform variable s < 0 satisfying 
s min < s < there is at least one r 0 at which the Laplace- 
transformed Lame parameters for Maxwell rheology. 


A(s, r) = 


AS + fixit] 
S A- pit] 


/<(s, r) = 


P s 

s + fi/rf ’ 


become singular. We will call the set of all r 0 a singular bound. 
For a model with N layers, there are only N points within the 
singular bound. We can prove (e.g. Fang, Hager & Herring, in 
preparation) that a mode can never be exactly at the singular 
value s= —fii*\ except when two adjacent layers degenerate 
into a single uniform layer. However, as the number of layers 
increases to infinity, the singular bound tends to span the 
entire closed interval [s^, s max ], and we have trouble 
determining if there are modes within the singular bound. This 
situation is equivalent to the 'theoretical crisis’ with the secular 
equation at N -> co. Following the spirit of using more and 
more layers to approximate a continuous function, we can 
interpret the singular bound [5^,5^^] as representing a 
continuous spectrum of modes (Han & Wahr 1995). Modes, 
as we call them, are indeed poles on the complex Laplace 
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transform plane. Poles must be isolated. A cluster point of an 
infinite number of poles is a non-isolated singularity (Colombo 
1983). As poles tend to be continuous, every point becomes a 
cluster point, and thus a non-isolated singularity. This is an 
unusual situation even in the context of complex analysis of 
one variable. It is not the intent of this paper to establish a 
generaMreatment in mathematical theory of such ‘continuous 
poles’. Rather, we will focus on how and to what extent the 
‘continuous poles’ contribute to the relaxation process of a 
Maxwell earth. This brings about the key question: how can 
we investigate the ‘continuous poles’? 

There are two major obstacles associated with the normal 
mode method that have previously prevented us from making 
a detailed investigation of the "continuous poles’. First, the 
normal mode method is designed on the assumption of isolated 
poles. In order to trap a mode, we need a small interval 
Oi, s 2 l If the determinants of the boundary matrix Mfo) and 
M(s 2 ) are of opposite signs, and M(s) is approximately a linear 
function in [s t , s 2 ], then there is a pole at s 0 (s t < s 0 < s 2 ) such 
that M (s'q ) — 0. This small interval is not only practically 
important, but also theoretically required as, again, a pole 
must be isolated such that there is a neighbourhood of s 0 in 
which no other poles exist. When poles become continuous, 
the distance between each adjacent pole tends to zero, and the 
basic assumption of isolated poles no longer applies. Even if 
we can find a way to circumvent the difficulty and manage to 
isolate the poles, the obstacle is still not removed. As discussed 
later on, since the number of poles is infinite, the residues of 
the poles must be zero, except for a finite number of them. 
However, ‘continuous poles’ as a whole may well have a non- 
zero contribution and may even form the entire contribution 
to the relaxation process (see Section 4 below). This problem 
reflects a fundamental difference between a discrete treatment 
and a continuous treatment. 

Runge-Kutta propagation is a method of using discretized 
layers to approximate a continuous function. One may reason- 
ably argue that we are practically dealing with discrete modes 
by using Runge-Kutta propagation while we are thinking of 
'continuous modes’. This argument is only one side of the 
story. On the other hand, discretization is a process of sampling. 
By sampling, we lose resolution power. If the original structure 
has two modes quite close to one another, a discretization 
with poor resolution power will lump them together to form 
a 'continuous spectrum’. This is quite similar to the Nyquist 
sampling problem in the area of data processing. What is 
happening on real numerical groups is a trade-off between 
these two conflicting mechanisms. We can look at the problem 
in a slightly different way. A discretization introduces error to 
the solution of a continuous structure. It is the same source of 
error as for a layered structure. A good numerical algorithm 
should be able to identify the differences between a continuous 
structure and a layered structure above the error level. Below 
the error level, it is meaningless to talk about either a continu- 
ous structure or a layered structure. We assume, throughout 
this paper, that our numerical algorithm is ‘good’. If this 
assumption were to be invalidated by actual calculations, we 
have to consider replacing the currently used Runge-Kutta 
propagation with some more accurate scheme — for instance, 
upgrading the Runge-Kutta scheme from fourth order to 
sixth order. 

The second problem with the normal mode method is that 
it forces the propagation to cross the singularity radius r 0 . The 


normal mode method used is an extension of the seismic 
normal mode technique (e.g. Takeuchi & Saito 1972). Starting 
from the centre of the Earth, numerical integration of the 
dynamic equations propagates upwards to the surface to meet 
the boundary conditions for a trial mode. There is never a 
problem with such propagation for an elastic earth model, but 
there is a problem for a Maxwell viscoelastic model in the 
Laplace transform domain; namely, the propagation breaks 
up at r 0 for a trial mode s within the singular bound. Forcing 
the numerical integration to go through r 0 results in a series 
of noisy singular solutions for the set of tnal modes s within 
the singular bound (e.g. Wu & Peltier 1982). These noisy 
singular solutions have caused concern (e.g. Mitrovica & 
Peltier 1992; Fang & Hager 1994; Han & Wahr 1995). Han & 
Wahr (1995) introduced a ‘zero crossing criterion’ to remedy 
this situation. It helps to identify some prominent modes that 
would otherwise have been buried among the noisy singular 
solutions, but it does not avoid the appearance of such noisy 
singular solutions. Thus it is of little help in assessing the 
singular bound contributions. This problem is related to 
the first obstacle mentioned above, namely, that the power of 
a discrete pole method becomes diminished when dealing with 
‘continuous poles’. 

For most practical problems, the Laplace transform is used 
to improve the convergence of transform integrals. For a 
Maxwell rheology, however, the convergence is not a problem 
at all, even for the ordinary Fourier transform. Fang & Hager 
(1994) realized that the singularities can be avoided by 
extending the ordinary Fourier transform into the complex 
frequency domain (equivalent to the Laplace transform in its 
original form) but the time inversion along this line is extremely 
difficult. Fang & Hager (1994) modified the inversion pro- 
cedure, and called it the complex real Fourier transform 
(CRFT) method. A more straightforward method that avoids 
the singularities is proposed by Hanyk et al, (1995). This is a 
direct time-domain solution so that singularities appear- 
ing in the frequency domain no longer exist. Yet none of these 
‘lumped response’ methods is able to provide direct 
answers concerning the roles that the singularities play in 
the relaxation process. 

In the past 20 yr, extensive investigations have been earned 
out on the discrete modes caused by the seismicallv identified 
discontinuities such as the 670 km boundary, CMB, the litho- 
sphere and so on. Luckily enough, for the simple layered 
viscosity models used in these investigations (e.g. Nakada & 
Lambeck 1989; Tushingham & Peltier 1991; Han & Wahr 
1995), these prominent modes caused by the discontinuities 
are mainly located outside the singular bound. Furthermore, 
these prominent modes seem to dominate the time responses. 
This latter conclusion comes from the ‘inviscid fluid criterion’ 
proposed by Wu & Peltier (1982). The ‘inviscid fluid criterion’ 
concerns the final state of equilibrium with time 
Consider two admissible respones w(f) and vv L (r ): 

w 1 (r) = w(t) + e" <r ' £ r n (ff)t n , cr>0. 

n = 1 

The a term in vv L (r) is a typical component for a second- or 
higher-order mode. The final states w(cc) and w^oc) are the 
same. In this sense, the inviscid fluid criterion is not unique 
and incomplete. Nevertheless, it does indicate that, for 
these specific layered viscosity models, the singular bound 
contribution does not extend to a sufficiently long time. 
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Studies of the creep mechanisms for silicates (e.g. Weertman 
1978; Ranalli 1991) indicate that a ‘realistic’ viscosity structure 
is likely to be piecewise continuous, and a layered structure is 
‘unrealistic’. Could the general features revealed from studies 
of those ‘unrealistic’ models apply to the ‘realistic models’? 
Could a continuous structure generate prominent discrete 
modes? Is there a possibility that a singular bound contribution 
dominates the time response? All these important questions 
among others are related to one, namely, what is the real 
picture behind the noisy singular solutions? 

In this paper we try to give a thorough and unified account 
of the singularity contributions to the relaxation process. This 
goal is achieved by the introduction of what we call the 
‘isolation function' in the Laplace s plane. Using the isolation 
function, we are able to isolate the effect of each individual 
singularity, whether it is a simple pole, part of a ‘continuous 
mode’ spectrum, a branch point, or even essential singularities, 
etc. Because our treatment is fairly general, and no presumption 
is made of the nature of the singularities, we will follow the 
informal spirit to use the terms ‘continuous poles’ or ‘continu- 
ous modes'. We start in Section 2 with the basic mathematical 
framework and an outline of some exact solutions needed for 
later analysis. Section 3 is devoted to the introduction of the 
isolation function. Then, in Section 4, we present the results of 
isolation functions for a number of admissible viscosity models. 
Also in this section, we investigate the contributions of the 
singular bounds to the total load Love numbers. Finally in 
Section 5, there is the conclusion. In the Appendix we give a 
procedure for dealing with the singularities in the search for 
discrete modes on the negative side of the real axis in the s 
plane. The key to the procedure is a power series solution of 
the basic equations expanded about a singularity point. We 
provide details of the solution. 


2 EXACT MODE SOLUTIONS 
2.1 Basic equations 

We consider, throughout this paper, a non-rotating, spherically 
symmetric, self-gravitating, incompressible mantle with uni- 
form density p m and uniform shear modulus u. The effect of 
an inviscid core is taken into account by the boundary 
conditions and density contrast Ap CMB at the core-mantle 
boundary (CMB). The only reason for choosing such a simpli- 
fied earth model is that we have a closed analytical mode 
solution for layered rheological parameters, so that a direct 
check can be made for the main results obtained numerically. 
For an infinitesimal displacement field u, we have the unified 
quasi-static governing equations in the Laplace transform s 
domain (e.g. Wolf 1994; Fang et ai in preparation): 

dE ( du \ 

- Vpi — p m V<t>i 4- — 2 — -f f x V x u I — £V x V x u = 0, 

(1) 

V • u = 0, (2) 

V 2 ^ = 0, (3) 

where p t is the incremental reduced pressure, <j> { is the 
incremental self-gravitational potential caused by the surface 
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and CMB defiections, and 


rp elastic solid , 

«s viscous tlu id , 

— . — viscoelastic bodv . 

Vp/r/ 4- s 

Following the conventional treatment (e.g. Alterman, Jarosh 
& Pekens 1959), we decompose u, and into harmonic 
expansions: 


«= z 


f U n [r, s)P n { cos 0) 4 6 — K,(r. s)P„(cos 0) 
ctl 


CO 

<t>i= Z <I U r - s)P„(cos 0), 

n - 1 

Pi= Z s >^( cos 0)’ 

n = 1 

where P n represent the Legendre polynomials, and 0 is the 
colatitude. The solution for (3) is straightforward; 

® n = A n (s)r" + B n ( S )r-''-\ 


where A n and B n are the integral constants to be determined. 
To solve for the rest of the unknowns, let us define the column 
vector 


Y = ([/„, t T„ - K, T e „) r 


(4) 


where T rn and T 0n are the normal and shear stress components, 
respectively, for harmonic degree n . Then, eqs ( 1 ) and (2) can 
be reduced to a set of ordinary differential equations for each 
harmonic degree n: 


d Y 

— - MY, 
dr 

where 


r 

12 E 
r 

1 

r 

6 £ 
r 2 


M : 


n(n - 1) 


6En{n 4- 1 ) 
r 
i 
r 

2£(2m 2 4- 2ji — 


nin 4- l) 
r 
1 

£ 

3 

r " 


(5) 


For those £ of complex value (e.g. Section 3), eq. (5) could be 
made double size: 


d 

yc 


m r 

-M,- 

~y r - 

Jr 

_Y,_ 



m r . 

. Y, _ 


where Y R , M R and Yj, M 1 are the real parts and imaginary 
parts of Y and M, respectively. 

Free-slip boundary conditions are imposed on both the 
surface and CMB (e.g. Richards & Hager 1984) for a pure 
relaxation process, while for the surface loading problem, 
additional loading terms introduced by Longman (1962) are 
added to the radial stress component of the surface boundary 
condition and the perturbed gravitational potential. Physical 
variables Y, <t>„, and d<t> n /dr are all continuous at viscosity 
discontinuities. The extension of the boundary conditions from 
real to complex is straightforward. 
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The viscous parts of the vertical and horizontal load Love 
numbers, h n and J n , for a Maxwell rheology are defined 
following Peltier ( 1974) as 

K = K-K> l = (7) 

where h n , /„ are the total load Love numbers and h c n , l c „ are 
their elastic counterparts. The third Love number, k n1 is strictly 
proportional to h n for our earth model, and is therefore not 
considered. For simplicity, we also call h„ and /„ the Love 
numbers. 


2.2 Two-layer solution 

We outline and discuss the results. Details of how these results 
are obtained can be found in Fang et ai (in preparation). 
Consider a set of two-layer viscosity models with a uniform 
lithosphere of 120 km thickness and increased viscosity > 7 Uth 
overlying a uniform mantle of a fixed viscosity j/ m at the value 
of the so-called background viscosity, r] 0 = 10 21 Pa s. Other 
necessary physical parameters are listed in Fig. 7. The time 
response of the Maxwell Love numbers to an impulse can be 
written as a sum of modes: 




exp(S;t), 


Sj < 0 , 


( 8 ) 


where Sj is the eigenvalue of the jth mode. Two out of the four 
modes are caused by the density contrasts at the surface and 
CMB, and the other two are excited by the viscosity discont- 
inuity at 120 km depth. Note that there are two discrete 
singularities for such tw r o-layer viscosity models, namely 


1 7lith ^ _ *Tm ^ g | 

'^max M ^min M 

It is more convenient in this situation to use the notion of 
reciprocal singular bound defined as [— l /5 tnui , — l/s ma x]> since 
— 1/s corresponds to the relaxation time. Fig. 1 shows relax- 
ation time spectra for six viscosity models for both a viscous 
and a Maxwell rheology. We adopt the notation from Wu & 
Peltier { 1982) of mantle modes (M) and core modes (C), while 
we call the third spectrum of modes, to which they gave the 
name lithosphere modes (L), the viscosity modes (V), because 
the honzontal eigenfunctions of these modes clearly show the 
viscosity contrast (Fang et ai in preparation). Our identifi- 
cation of modes is also different from Wu & Peltier (1982) in 
the C modes and V modes before ‘transition degree’ 6. 

There is a fourth spectrum of modes which we name the 
surface modes (S). These modes have been practically neglected 
in previous studies but turn out to be very useful in our 
analysis. In fact, a Newtonian viscous fluid always has two 
groups of gravitational modes: the core modes (C) excited by 
the CMB, and the surface modes (S) generated by the surface 
boundary. For a uniform viscosity model, a Maxwell body 
also has only two spectra of modes, S and C (Fig. la). The 
overlapping M and V modes in Fig. 1(a) are zero amplitude 
modes, as there is no viscosity contrast there. Non-zero M 
and V modes occur when the viscosity contrast is not zero. 
Fig. 2 shows the relative strengths of the C, V, and M modes 


in the Love numbers for the viscosity model in Fig. 1(f). The 
strength function is defined by 


kYj ' 


~h r V S j 
I l(t r >/ s ;l 

/= i 


i Yj — 4 


— i r" / s ■ 
i r 


j= 1,2, 3,4. 


I i.^i 

i=i 


This strength function is similar to that used by Wu & 
Peltier (1982). In their case, the signs of all trie amplitudes of 
the vertical responses are consistent, so that there is no need 
to take special care of the denominator. Here we encounter a 
situation where the horizontal responses can have opposite 
signs (see below). A direct summation will result in cancellation 
of the magnitude at the denominator. We choose the absolute 
value of the amplitude of each mode to form the denominator. 
The meaning of the strength function in this case is different 
from that used by Wu & Peltier (1982). 

As seen in Fig. 2, for models with a lithosphere thickness of 
120 km, the total strength of a Love number is mainly shared 
by the three M, C, and V modes. The contribution from the S 
modes is negligible. This is why the S modes have always been 
ignored in previous studies. However, this is not always true 
for different lithospheric models. If we (artificially) increase the 
thickness of the lithosphere to 1000 km, the picture of the 
strength function is quite different (see Fig. 3). 

Two important observations can be made concerning Figs 
L 2 and 3. First, the eigenvalues Sj of the two modes excited 
by the viscosity contrast are very close to each other when the 
viscosity contrast is small. It is reasonable to extrapolate that, 
as s min"^max> the g a P between the two modes excited by the 
viscosity contrast will tend to zero. Additional numerical tests 
support this extrapolation and we do not see how it could 
happen otherwise. This observation implies both the existence 
of ‘continuous modes’ for a continuous viscosity and the 
association of the singular bound with the ‘continuous modes’. 
Note that a viscosity contrast corresponds to two discrete 
singularities (see eq. 9), and a viscosity contrast creates two 
modes. Therefore, we can say that one singularity claims one 
mode, or equivalently, the number of the V modes and M 
modes is equal to the number of singularities. This singulanty- 
mode-correspondence feature extends to layered models of any 
number of layers (Fang et ai in preparation). 

Secondly, if we fix the thickness of the lithosphere at 120 km 
and increase the viscosity contrast, or, equivalently, increase 
the reciprocal singular bound, more modes will be trapped 
into the reciprocal singular bound. Fig. 2 shows that, for 
sufficiently high viscosity contrasts (larger than a factor of 20), 
the strengths of the Love numbers mainly come from the 
modes trapped within the reciprocal singular bound. On the 
other hand, if we fix the reciprocal singular bound and increase 
the thickness of the lithosphere, the strength of the Love 
numbers will quickly leak to the modes originally outside the 
reciprocal singular bound (Figs If and 3). 

These two observations provide a basic and intuitive guide 
to the rest of the analysis. 


3 THE ISOLATION FUNCTION 
3.1 Definition 

Consider an open-ended box contour C in the complex z = 
x 4- iy plane (Fig. 4). Fixing ail the dimensions of C except for 
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10 ° 10 1 10 2 





Harmonic degree 


Figure 1. Relaxation tjmes of the two-layer viscosity, otherwise uniform, mantle. The thickness of the lithosphere is fixed at 120 km, the mantle 
viscosity rj m is fixed at 10 21 Pas, and the lithospheric viscosity ^ lith varies as indicated in each of the windows. Other necessary physical parameters 
needed for the calculations are listed in Fig. 7. The circles are for a Newtonian viscous mantle, and the dots are for a Maxwell viscoelastic mantle. 
The reciprocal singular bounds are marked by dashed lines (from Fang et ai in preparation). 


the x coordinate of the open end, one can create a function of 
x by means of the contour C integral of a specific integrand 
w ( z ) . Next, we shrink the height e of C, which is parallel to the 
imaginary axis, to zero in such a way that the real axis is 
always between the two sides of C ( Fig. 4). We detine the 
isolation function / w (x) as 

IJx) = iim — - (b w(z) dz , (x < .x fix ) . 

£-*0 7U J C(x) 

Except when mentioned otherwise, we assume that singularities 
of w{z) are on the real axis and that poles are all of first order. 


One can verify that, if w(z) has one pole at x 0 within the 
interval ( x , x fix ), the isolation function I w (x) is equal to the 
residue of the pole. In fact, as the height shrinks to zero, the 
contribution of the height on either side of a closed box 
contour to the integral will tend to zero; as a resuit, the closed 
box is equivalent to an open-ended box in the limit. As a 
demonstration, we examine the isolation function of the 
integrand 

Mz)= — - — , 

z — z 0 
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Figure 2. Percentage strength of the modes in the viscous part of the 
load Love numbers. The definition of the strength can be found in the 
text. The earth model use in the calculation is in Fig. Iff). The S mode 
contribution is not plotted. The total contribution of the plotted C, V 
and M modes is nearly 1 CX> per cent of the total strengths of the Love 
numbers, indicating that the S mode contribution is negligible. 


where z 0 = x 0 , (,x 0 < x fix ). One can easily obtain 


/„(x) = lim — 


l,m-U 

'~o 2m J c 


. 1 
dz = lim - 

E — 0 Jl 



8 

(X - X 0 ) 2 4- £“ 


dx' 


x < x 0 , 
.x>.x 0 . 


.x fix >0. 


( 10 ) 


This is a step function, jumping in value at the pole z 0 of w. If 
there are several poles on the left side of ,x fU , I„(x) will be a 
muitistep function with the same number of jumps (falls) as of 
poles (Fig. 4). The amplitude of each relative jump (fall) is the 
residue of the pole. If there is a spectrum of ‘continuous poles’, 
or a branch cut, I w (x) will become a continuous function of x, 
as illustrated in Fig. 4. Note that eq. ( 10) indicates that the 
limit sign and integral sign in the isolation function are 
generally not interchangeable. 


3.2 Continuous poles 

The mechanism for causing the continuity of the isolation 
functions by branch cuts is well known (e.g. McLachlan 1955). 
But ’continuous poles* or a continuous distribution of non- 
isolated poles are not often seen in the literature. Thus, we 
need to formulate them here. Following an informal spirit, a 
continuous viscosity can be viewed as the limit of a sequence 
of multilayered viscosity models, with the number of layers 


Figure 3. As Fig. 2, except that the thickness of the lithosphere is 
increased to 1000 km. The total contributions of the plotted C, V and 
M modes are significantly less than 100 per cent of the total strength 
of the Love numbers, especially at higher degrees. The deficiency of 
the total strengths is the contribution of the S modes. 
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Figure 4. A picture demonstrating the open-ended box contour and 
the isolation function, e is the height of the contour and is supposed 
to shrink to zero. No singularities, except the four poles represented 
by the four dots, are on the real axis between .x and .x rix . The dashed 
line indicates that each pole corresponds to a step in the isolation 
function I(x). The continuous curve illustrates what the isolation 
function looks like when branch cuts, or ‘continuous poles’, occur 
between .x and .x fU . 
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tending to infinity. Without seeking mathematical generality, 
we can formulate ‘continuous poles’ as the limit of a sequence 
of finite poles. 

Suppose there are M poles of the integrand w(z) = /(z)e rf 
within a bounded interval (x 0 ,x fix ). The isolation function 
with a constant parameter t is 

U.x, f ) = I r(xj)e x i' !f(Xj - x ) (x 0 < x, < ■ < x M < x fi J , 

y=t 

wTere r(xj) is the residue of the pole at x } , and H is the 
Heaviside function. Since (x 0 , x rix ) is bounded, we will encounter 
cluster points in the sequence of poles as M~*o o. If the number 
of cluster points is finite, we have a dense mode distribution 
only at these isolated cluster points. This is a peculiar situation 
where we find neither a discrete mode distribution nor a 
continuous spectrum of modes. To form a continuous spectrum 
of modes, the number of cluster points must be infinite. For 
simplicity, we assume that the poles are always nearly evenly 
distributed within (,x 0 , x fix ). This assumption is supported by 
the observations from the experiment with simple layered models 
(see Section 2). For a distribution with an infinite number of 
poles in a finite region, the residues of all the poles must tend 
to zero as M — *■ go, except for a finite number of them, otherwise 
0 will tend to infinity, which is physically unacceptable. In 
mathematical terms, we have that 

lim r{xj) = 0, lim (x ) + 1 -x.) = 0, ;=1,2, 3,... 

■V— co M oo 

(11) 

are true almost everywhere within (x 0 , x fix ). 

Using eq. ( 1 1 ), we can define the density of the residue <j( x i ) 
as 


Z(Xj) ■ 


r(xj) 


(X }+l -Xj)' 

Note that, at those exceptional points Xj where the residues 
have finite amplitude lim w _^r(x,) ^0, the density s "(x ; ) is in 
the form of a Dirac delta function. So, generally, we have 

M 

Ux,t)= lim V i(xj)e*>' H(x, — x)(x i+ 1 — X,) 


cOOe* f dx ' . 


(12) 


Relations (11) and (12) give the meaning to the frequently 
used term ‘continuous spectrum of modes’ in this specific 
situation. In particular, relations (11) simply state that the 
residue of almost every individual pole in the continuous 
spectrum is zero, while eq, ( 12) states that the zero-contnbution 
poles as a whole may have a significant contribution to the 
isolation functions. The treatment introduced here is not 
limited to the first-order poles going continuous. The difference 
between first-order and higher-order poles only lies in the 
method of extracting the residues r(xj). For second- and higher- 
order poles, the residues r are also, in general, functions of 
time. In our formulation for the continuous spectrum of modes, 
the time variable can be considered as a constant. For a 
distribution of non-pole singularities, the formulation of a 
continuous spectrum of modes from ‘continuous poles’ may 
fail, but the isolation function as defined by a contour integral 
on the complex plane is still meaningful. In fact, we can use 
the isolation function, at least symbolically, to define a continu- 
ous spectrum of modes (see below). 


3.3 The isolation function for Love numbers 


The endpoints of the singular bound [ 5 min ,s max ] of a 
continuous viscosity model are 


/' * 


/‘ 

n(r)_ 

| - = max 

. r r( r ).. 


For a continuous viscosity, every single point within [s min , s max ] 
corresponds to at least one singularity while for a two-layer 
viscosity, only the endpoints of the singular bound, s min and 
correspond to singularities. However, distinguishing 
between the two types of singular bounds does not seem to be 
important, and we will not do so. 

The Laplace inversion for the Love numbers in response to 
a Heaviside load history is 


vv n (r)= lim — 
.v-oo 2 m 


'*<r +■ * jV 
Ja — iS 



ds, 


a > 0, 


(13) 


where \v„ represents either h„ or l n . Performing the usual trick 
in complex analysis, the integral (13) can be carried out, in 
the s plane, on the open-ended box contour (Fig. 5). Note that 
the contour in Fig. 5 excludes the origin s = 0. The right end 
of the contour extends a little bit to the right of s max bv a 
small distance S l >0. We shrink the height of the contour the 
way we did before and reduce the complex expression (13) to 
a real expression 


1 

w n {t) = lim- 

71 



dco 


+ X-(e^- 1), 

J S J 


(5 = — 03 -r is) , 


fJ.S) = 


\vU-a) sin et - e cos gr) 4- cos et + e sin et) 

e 2 + or 


(14) 


Im 



Figure 5, The contour for the Laplace inversion. The radius of the 
arcs will tend to infinity, and the height of the open-ended box shnnks 
to zero. The isolation function contour is shown as part of the 
shrunken open-ended box (the thick segment}. As a comparison, the 
integral path used in the CRFT method (Fang & Hager 1994) is 
also indicated. 
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where vvj, and vv‘ n are ihe real part and imaginary parts of vv„($), 
respectively, the integral limits are oi min = — s max — S lt co max = 
-s mm + d 2 , d ly S 2 > 0 and r Jf Sj are, respectively, the residual 
and location of the yth pole found outside the integral interval 
-ojmin). In reaching eq. (14), we have used the prop- 
erty that the complex conjugate of the Laplace transform of a 
real function can be obtained by replacing the Laplace variable 
s with its conjugate. Lq. ( 14) can be interpreted as the complete 
eigen-expansion for the Love numbers with both a continuous 
spectrum and a discrete spectrum. 

Again, the formulation of the integral term in ( 14) does not 
exclude the possibility of a number of modes with finite 
amplitudes sticking out among other 'continuous modes’. In 
this case, the density function f n (s) will tend to the Dirac delta 
function as e-*0. Another possibility is that there are branch 
cuts extending out of the singular bound. In this case we can 
use the quantities S { or S 2 , or both, in case it is necessary, to 
extend the singular bound so as to include the branch cuts in 
the integral ( 14). If no branch cuts extend out of the singular 
bound, we can set and S 2 to small positive values to make 
sure that the endpoints of the singular bound are fully counted 
in a numerical evaluation of the integral (14). The variable 
co represents the reciprocal relaxation time. It also conforms 
to the general sense of a function with a positive ascend- 
ing argument. For simplicity, both (— cv max , — co mxn ) and 
(cu rain , to ma J are referred to the same extending singular bound, 
as they indeed are. 

The isolation functions for the Love numbers are defined 
on the extended singular bound. With a slight change of 
notation, we have 

1 

f* ( — co, f} = lim - 

(cu mm < tt> < w m , x ) . (15) 

One question remains: how does the isolation function 
behave numerically? The key issue is whether a numerical 
integration of (15) can give a robust approximation to a step 
function. Fig. 6 shows the numerical result of the test integral 
( 10). Since the numerical integration will blow up at the exact 
value e = 0, the limit of £ is stopped at z — 0.002. The integral 
spacing is chosen as 0.001, and the upper limit x fix is fixed 
at 1. Fig. 6 clearly indicates that the isolation function is robust 
enough numerically to isolate the effects of modes. 


f„{s')e ait <k o'. 
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Figure 6. The numerical result of the isolation function for the 
integrand 1 jz. The uniform spacing for the numerical integration is 
chosen as 0.001. and the limits stopped at e = 0.002. 


4 RESULTS 


4.1 Viscosity models 

Four sets of viscosity models. A, B, C, and D, are used in this 
study (Fig. 7). The singular bounds s max ] for all of the 

viscosity models here are the same, with the endpoints 


_iL . F 

^mirv » "'max r~> ' 

no 53.7a; 0 

Obviously, the singular bound is not a sufficient measure to 
count for the contributions of various viscosity structures. We 
define the convexity of the viscosity models as illustrated by 
Fig. 7. The exact measure of convexity for a viscosity model 
n(r) is 


sup 

r CMB ^ r ^ r 610 


nCMB n 670 ( \ t \ 

^670 + : — (r - ' 670 ) - n(r) 


for A , 


sup 


n(r) - ^70 + 


r CMB r 670 

n a - n 670 


r a - r. 


(r — r 670 ) 


(16) 


for B, 


where subscript a denotes the surface and other subscripts 
mark the depth. The definition (16) simply states that, within 
the same singular bound, a linear variation of viscosity always 
has zero convexity, the quadratic and exponential variation of 
viscosities always have negative convexities and a layered 
viscosity always has positive convexity. Note that the convex- 
ity for a layered viscosity is just the thickness ol the layer 
(see Fig. 7). 


4.2 Isolation functions for different viscosity models 

The main results of this paper are the isolation functions (14) 
for the viscosity models in Fig. 7. We pick up harmonic degrees 
4 and 15 for presentation here. The times are fixed at 1.2 kyr 
for degrees 4 and 4.2 kyr for degree 15. The horizontal isolation 
functions are in the form of 

nl\ n (—(D, t ) . 

Fig. 8 shows the isolation functions for lower mantle 
models A. Surprisingly, all the models in A have discrete mode 
structures within the singular bound. With the aid of Figs L, 
2, and 3, we can easily identify the M, V and C modes. The 
inclusion of the lithosphere brings about an additional viscosity 
discontinuity at the bottom of the lithosphere which will 
maintain two more modes. In total, the layered model in A 
has six modes. Because it is not a goal of this paper to identify 
modes, we simply term the modes other than M, V and C as 
X t , X 2 , and X 3 , etc. For the continuous viscosity models in A, 
the minor modes are smoothed out, but the major modes 
remain. An inspection of the curves in Fig. 8 shows that the 
remaining major modes are also slightly smoothed, suggesting 
that the effect of the continuity of the lower mantle viscosity 
is not strong enough to split the discrete mode contribution 
into a continuous mode distribution. On the contrary, the 
continuous structure of the upper mantle viscosity has much 
stronger effects on splitting the discrete mode contnbution 
into a continuous mode distribution. This feature can be seen 
in Figs 9 and 11. Furthermore, Figs 9, 10 and 11 quantitatively 
illustrate the situations where no discrete modes exist within 
the singular bound. 

There have been a number of studies, some very recent, on 
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Viscosity models 


Figure 7. Viscosity models. Both model sets A and B contain a two-layer, an exponential, a quadratic and a linear model, while C and D do not 
have the two-layer model. Viscosity variations of A are within the lower mantle, starting with the value r\ m = at 670 km depth to the value r} m = 
53. 5^0 at the CMB. The upper mantle viscosity of A is uniform at the background viscosity rj 0 except for a lithosphere. The lithosphere of A is 
modelled as a 120km thick layer with the viscosity exponentially increasing from ij 0 at the bottom to 10 6 >j o on the surface (the shaded area). 
Viscosity models in B do not have a lithsophere and their variations are within the upper mantle, starting with the value - rf 0 at the 670 km 
depth and increasing to - 53 at the surface. The lower mantle viscosity of B is uniform at the background viscosity rj 0 . The model set C is 
a combination of A and B without the lithosphere, while D is a combination of A and B with the lithosphere. 


the modal branches associated with simple layered viscosity 
models (e.g. Wu & Peiter 1982; Wolf 1985; James 1991; 
Amelung & Wolf 1994; Han Sc Wahr 1995; Mitrovica & Davis 
1995). In all of these previous analyses, the MO mode (equival- 
ent to the M mode in this paper) is termed the fundamental 
branch of relaxation, and interpreted as arising from the 
density discontinuity at the outer surface of the Earth. This 
physical description of the MO mode hardly fits the results of 
our analysis. Because the modal branches exchange elements, 
at lower degrees (n < 6), it is hard to decide, according to the 
‘movie’ in Fig. 1, whether the MO mode is related to the 
surface density contrast or to the viscosity contrast at 
the bottom of the lithosphere. At higher degrees (n > 40), it is 
clear, by comparing the viscous mode branches and viscoelastic 
branches in Figs 1(e) and (f), that the M0 and L0 modes 
(equivalent to the M and V modes, respectively, in this paper), 
belong to the pair of modes arising from the viscosity disconti- 
nuity at the bottom ol the lithosphere. The mode arising from 
the surface density contrast joins the S branch at degrees 
n > 40. When the viscosity contrast at the bottom of the 
lithosphere is very large ( Figs le and f), the S modes are way 
above the rest of the modal branches, and the contribution of 
the S modes to the time response is negligible (Fig. 2). 
Furthermore, if the M0 mode constitutes fundamental modes 
arising from the surface density discontinuity, it should stand 
out at all the harmonic degrees and display a certain kind of 
independence from the viscosity structure. Fig. 9 merely shows 
that the M mode at very low degree (e.g. degree 4) does show 
some independence from the viscosity models, suggesting that 
the M mode is due to the outer surface density contrast. 
However, at degree 15, the M mode is completely split into 
continuous mode distributions by stretching the viscosity 
contrast into continuous viscosity structures ( Fig. 9). This 
feature hardly supports the interpretation that the M0 modes 
are due to the surface density contrast at degrees n>15. 
Rather, it supports the assertion made above, that the M0 
modes and L0 modes at higher degrees (n > 15) are the twin 


sisters arising from the viscosity contrast at the bottom of the 
lithosphere. The validity of this assertion may even apply to 
lower harmonic degrees. 

A close inspection of the linear curve in Fig. 9 shows a 
number of blunt saw teeth in the vicinity of where the M mode 
arising from a layered structure is found. These saw teeth can 
be used to demonstrate how the M mode is split into a 
continuous distribution. A layered structure which has a 
positive convexity (the layered model in B in Fig. 7) generates 
the M mode. As we stretch the viscosity contrast into a 
continuous structure, the convexity of the model decreases. At 
zero convexity, we observe saw teeth in the vicinity of the 
former M mode, indicating that the single M mode is split 
into a number of small modes. As the convexity of the model 
decreases further, becoming negative, the saw teeth disappear, 
and the M mode is completely smoothed out. This phenom- 
enon implies that, for certain classes of viscosity structures (e.g. 
monotonic functions within the same singular bound), the 
transition from a discrete mode distribution to a continuous 
mode distribution could be gradual and smooth. 

For viscosity models having non-monotonic variation with 
depth, the situation is more complicated. A specific value in 
the singular bound cu(cu min < a> < oi max ) corresponds to at least 
two singularities. For the combined models in C and D we 
have one singularity in the lower mantle and another in the 
upper mantle. The combined effect of the upper mantle continu- 
ous viscosity and lower mantle continuous viscosity provides 
smoother (Figs 10 and 11) isolation functions, and no saw- 
tooth structure is found with the viscosity of linear variation 
with depth. The smoothing effect is largely due to the upper 
mantle continuous viscosity structure (see analysis of Figs 8 
and 9 above). The effect of the lower mantle continuous 
viscosity can be seen by comparing Figs 9 and 10. The pattern 
of gradual evolution in the isolation function from the single 
M mode structure to continuous mode distributions with a 
decrease of convexity disappears in Fig. 10, and a kind of 
irregular relationship in terms of their convexities shows up. 
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Reciprocal relaxation time co (1/ky) 


Figure 8. The isolation functions of the viscous part of the vertical 
and horizontal load Love numbers for viscosity models in A at 
harmonic degrees 4 and 15, respectively. To demonstrate the mode 
structures better, each curve is shifted consecutively by the amount of 
the annotation interval. Since the isolation functions all quickly tend 
to zero, one can easily identify the shift of each curve and calibrate it 
at the right zero level. 

In general, all these figures (Figs 8 and 11) display a compli- 
cated relationship between a viscosity model and the spectra 
of the relaxation eigenfunctions associated with the model. 
This is in sharp contrast to the relationship between an elastic 






Figure 9. The same isolation functions as in Fig. 8, but for viscosity 
models in family B. The plots here are not shifted. 


earth model and the elastic eigenfunctions (normal modes). 
The distribution of eigenfunctions for all elastic earth models 
is the same: grand modes 0 So>oSi an< ^ overtones (e.g. 
Lapwood & Usami 1981). The eigenspace for all admissible 
elastic models can thus be easily established. However, here 
the distribution of relaxation eigenfunctions associated with a 
viscosity model depends upon the viscosity model itself in a 
non-linear way. This situation makes it much harder to estab- 
lish the relaxation- eigenspace for all admissible viscosity 
models. 

The isolation functions (Figs 8 and 11) provide another 
means of revealing the differences between the vertical and 
horizontal responses. The difference can be observed with 
either the success (Fig. 8) or the failure of the discrete mode 
structure (Figs 9, 10 and 11). The essence of the difference, 
using the mode language, lies in the fact that major modes of 
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Figure 10. As Fig. 9, but for viscosity models in C. Figure 11. As Fig. 9, but for viscosity models in D. 


the vertical response all have the same sign, while for the 
horizontal response, the V and M modes have opposite signs 
(Figs 2 and 3; Fang & Hager 1994; also see Mitrovica & 
Davis 1995). 

4.3 Contributions from singular bounds 

We have just seen in the last subsection both discrete mode 
and continuous mode structures within the singular bound. 
The following questions arise: (i) how much does a singular 
bound contribute to the total strength of Love numbers; 
(ii) when can a singular bound be dismissed? We look into the 
problem by examining the time variation of the Love numbers. 

We calculate the time variation of the Love numbers, which, 
according to the inviscid fluid criterion (Wu & Peltier 1982), 
should be independent of viscosity models as f->oc. We then 
calculate the time vanation of the singular bound contributions 


to the Love numbers. A comparison between the two indepen- 
dent calculations will provide a direct assessment of how much 
a singular bound contributes to the Love numbers at a given 
time. Since the horizontal motion is not sensitive to gravity, it 
is not directly related to the mechanism of isostasy. For 
this reason, we save the discussion of the horizontal motion 
for elsewhere, and concentrate on the vertical motion here. 
Figs 12 to 15 show the results for vertical responses. The time 
variations of the Love numbers are calculated using the CRFT 
method (Fang & Hager 1994). The singular bound contri- 
butions are calculated with the isolation function ( 15) by fixing 
co at its lower limit, and letting the time f vary. 

A striking feature in Figs 12 to 15 is that the contributions 
of the singular bounds to the total Love numbers depend 
upon, and are even quite sensitive to, the convexities of the 
viscosity models (remember the singular bound for these 
viscosity models is the same). The characteristic relaxation 


0 1995 RAS, GJl 123, 849-865 




860 M. Fang and B.H. Hager 





Figure 12. Time variations of the total Love number and the singular bound contributions for the viscosity models in lamily A. The total Love 
numbers are calculated using the complex-real Fourier transform (CRFT) method (solid line). The singular bound contributions are calculated 
with the isolation function (dashed lines). 


time of a Maxwell rheology is which measures the effect 
of viscous damping (resistance) of the dashpot element against 
the instantaneous elastic response of the spring element. The 
convexity of a viscosity model as a whole measures the relative 
importance of the viscosity effect among the models within the 
same singular bound and thus represents the relative degree 
of viscous damping. With the singular bound fixed, the larger 
the convexity, the heavier the viscous damping. As revealed in 
Figs 12 to 15, it takes a longer relaxation time in the case of 
heavier viscous damping for the total Love numbers to reach 
the final state of complete isostasy. The contribution of the 
singular bound to the total Love numbers, on the other hand, 
is smaller in the case of heavier viscous damping. For the 
exponential viscosity models, which have the smallest con- 
vexity, the singular bound contribution dominates the strengths 
of the total Love numbers. We have seen in the last subsection 
that there is effectively no fundamental modal branch arising 
from the outer surface density contrast and spanning the 
spectrum of all harmonic degrees for all viscosity models. Here, 
we further confirm that, for a certain class of viscosity models, 
the contribution from the fundamental modes, if they exist, is 
negligible. 


The layered viscosity model in group B has an exact mode 
solution, as outlined in Section 2. It is equivalent to a uniform 
lithosphere of 670 km thickness overlying a uniform mantle. 
As seen in Fig. 13 (layered), the modes trapped within the 
singular bound, including the M modes, no longer dominate 
the total strength of the Love numbers with such an artificially 
thickened lithosphere. The same feature has been observed in 
Fig. 3. Fig. 13 clearly indicates that what we have observed in 
Fig. 3 is just a special case of the general situation where, with 
the singular bound fixed, heavier viscous damping results in 
less singular bound contnbution to the Love numbers. 

The relationship revealed in Figs 12 to 15 between the 
convexity of a viscosity model and the singular bound contri- 
bution to the Love numbers is very useful in evaluating the 
previous results and predicting the singularity effect on a new 
model. As an example, we see that the predicted effective 
viscosity, based on diffusion creep theory, is close to an 
exponential function (e.g. Ranalli 1991). For an estimate, we 
choose the viscosity near the D" region to be 150^ o , and the 
viscosity at 670 km depth to be r/ 0 . These figures are quite 
conservative according to previous studies (e.g. Ranalli 1991). 
If the diffusion creep law dominates the lower mantle, we will 
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have a "realistic’ viscoelastic parameter rj/j.i quite similar to 
that of the exponential model D in Fig. 7 [the 'realistic' shear 
modulus, e.g. Dziewonski & Anderson (1981), also has a 
modest increase with depth]. Hence, we can predict from 
Fig. 15 that a significant part of the total strength of the Love 
numbers will come from the singular bound. On the other 
hand, the viscoelastic parameters in Tushmgham & Peltier’s 
(1991) model have nearly zero convexity and a very small 
singular bound. We can reasonably anticipate that the contri- 
bution of the singular bound to the Love numbers for this 
model will be small. These two predictions are indeed correct 
(see Fang et al. in preparation). 

A rule of thumb for spherical harmonic expansions in a 
wide range of problems is that the variation of upper mantle 
physical properties has little effect on the lower harmonic 
degree, while the variation of the lower mantle physical proper- 
ties has little effect on the higher degrees. This rule of thumb 
also applies to our problem. As seen in Figs 12 and 13, the 
lower mantle viscosity variation has little impact on degree 15, 
while the upper mantle viscosity variation has little effect on 
degree 4. Interestingly, at degrees where the viscosity variations 
do not have much effect, the singular bound contributions 
become as large as even up to 95 per cent of the total strengths 
of the Love numbers. For the combined models C and D, in 


which the upper mantle viscosity and lower mantle viscosity 
vary simultaneously, all degrees are affected as expected (Figs 
14 and 15). The lithosphere is a very shallow structure and 
therefore the C and D groups are nearly identical at degree 4 
and show some difference at degree 15 (Figs 14 and 15). 

5 CONCLUSION 

We have designed and implemented an isolation function for 
the load Love numbers by manipulating the topology of the 
contours in the Laplace transform s plane. Using the isolation 
function, we are able to scan through the entire singularity 
distribution created by a radially inhomogeneous Maxwell 
viscoelastic structure, and to touch upon a number of unsolved 
problems associated with the singularities. We sum up, in the 
following, the major results of this investigation. 

The discrete mode assumption for the Maxwell viscoelastic 
rheology is valid with finite layer viscoelastic parameter models 
( VEPMs). For VEPMs that are continuous functions of radius, 
the notion of discrete modes is only valid outside the singular 
bound. The general form of the dynamic response of a Maxwell 
viscoelastic earth contains a continuous spectrum of relaxation 
within the singular bound and a number of possible discrete 
relaxation modes outside the singular bound (eq. 14). 
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The discrete mode-like response can remain within the 
singular bound provided that the upper mantle VEPM is 
homogeneous or layered or the lower mantle effect dominates 
the viscoelastic relaxation [see also Hanyk et al. (1995) for the 
latter case]. In these cases, a discrete mode approach can be 


adopted as an approximation. However, the search for modes 
within the singular bound has to follow the correct procedure 
given in the appendix of this paper. The continuity of the 
upper mantle VEPM plays a decisive role in splitting the 
discrete modes into continuous mode distributions within 
the singular bound. The inclusion of the lithosphere has little 
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impact against mode smearing by the continuous upper 
mantle VEPM. 

The contribution of a fixed singular bound to the total 
strengths of the load Love numbers depends upon, and appears 
quite sensitive to, the convexity of the VEPM, which physically 
represents the relative degree of viscous damping among 
models within the same singular bound. For VEPMs of very 
positive convexity or heavier viscous damping, the load Love 
numbers are mainly determined outside the singular bound. 
For VEPMs of very negative convexity or lighter viscous 
damping, the major contribution to the load Love numbers is 
from within the singular bound. In particular, the theoretically 
predicted viscosity models have profiles of negative convexity, 
for example exponential functions of radius. We can be certain 
that singular bounds are indispensable in the studies of 
'realistic’ VEPM. 

We have observed elsewhere (Fang & Hager 1994; Mitrovica 
& Davis 1995) that, for layered VEPMs, the major modes of 
the vertical response all have the same sign, but this is not 
true Tor the horizontal response. We found in this paper that 
such an observation extends to the 'continuous modes’ within 
the singular bound (Fig. 9, model B, Fig. 10, model C). An 
indirect conjecture about such extension has been made by 
Fang & Hager (1994). The isolation function provides a direct 
and quantitative demonstration for the first time. 

The modal branch M, often called ‘fundamental’, is not 
generated by the outer surface density contrast alone, as was 
widely believed in the past. At very low degrees, the M modes 
demonstrate characteristics of being generated by the outer 
surface density contrast. At higher degrees (n>15), the M 
modes appear to be one of the twin modes generated by the 
viscosity contrast at the bottom of the lithosphere. When the 
viscosity discontinuity at the bottom of the lithosphere 
stretches into a continuous structure, the M modes split into 
a continuous mode distribution at higher harmonic degrees. 
The complexity of the physical causes of the M modal branch 
(also other modal branches) can be explained by the 'movie 5 
in Fig. 1 which shows that modal branches undergo an 
exchange of elements during the evolution from low- to high- 
viscosity contrasts. 

Finally, we discuss the validity of the simplifications made 
in this study. There are no new mathematical features involved 
in the extension of the incompressible uniform earth model 
used here to a ‘realistic’ earth model (e.g. PREM of Dziewonski 
& Anderson 1981). The major conclusions of this work can 
be extended to a realistic’ earth model. The isolation function 
(15) is based on the Heaviside loading history. If we use a 
smoother loading history, the isolation function should be no 
more 'bumpy 5 than that of a Heaviside, Therefore, our con- 
clusions can be extended to a 'realistic’ loading history. Our 
results concerning the convexity of the VEPM do not include 
piecewise continuous structures and oscillating structures of 
the VEPM. In these cases, the definition of convexity of a 
VEPM profile becomes meaningless. However, we can evaluate 
the degree of viscous damping by looking at the average of 
the viscosity structure. We believe that what we have found 
here, namely the heavier the damping, the less the singular 
bound contribution to the Love numbers, is true in general. 
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APPENDIX 


The key to the correct procedure of mode searching within a 
singular bond is to provide initial values near the singularity 
point r 0 . This requires an approximate solution of the basic 
eq. ( 1) expanded about r 0 . Applying the V x operator to eq. { 1 ) 
and substituting the harmonic expansion of the variables into 
( 1) and (2) we obtain 
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where for Maxwell rheology the coefficient E can be written as 


The derivation of the coefficients a k , h k , c k , d k is very messy for 
large k. This proves to be the major problem for high-order 
approximation. However, it is not so difficult to derive the 
first few terms. The crucial terms are the zero-order coefficient 
ti 0 , b 0 , c 0 and d 0 . They can be evaluated by the following 
limits: 
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Substituting eqs (A7), (A8) and (A9) into { A6) and making use 
of the singularity point eq. (A4), we obtain 
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Since the solution of (Al) itself is a meaningful subject and 
has not been tned before, in this appendix we first derive the 
power series solution of (Al) expanded about r 0 . The correct 
procedure for mode searching within the singular bound will 
be outlined afterwards. Here we only consider a Maxwell 
rheology. The other two rheologies, Newtonian flow and 
elasticity, do not suffer from the non-zero singularity point r 0 , 
and can be treated in a more complete fashion than we do 
here (e.g. Wu & Yuen 1991). 

Consider the solution of (Al ) within the mantle r CMB < r < r a . 
The only possible singularities of (i = 0, l, 2, 3) in the mantle 
region are at those r 0 values satisfying 

p. 4“ STfir^) = 0 , (^cmb — — r a) * (A4) 

It is easy to prove from eqs (A2), (A3) and (A4) that the 
following functions, (r — r 0 ) 4 ~ k q k (7c = U 2, 3, 4), are all regular 
at r 0 . Hence, we have the Taylor expansions 
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The general solution of (16) expanded at r 0 can be written 
as 

U n (r,s) = [_r-r 0 (s)y £ w k(s)lr “ r o(s)]*> (All) 

k — 0 

where X is the index for avoiding branch point. The parameter 
s is kept constant in seeking U n as a function of r. It will be 
dropped for simplicity. Substituting eqs (A6) and (All) into 
(Al), which is already multiplied by (r — r Q ) 4 , we obtain the 
indicial equation and recursive relation for the unknown 
coefficients w k : 

X(X — 1 )(X — 2)(a -3)4- d 0 X(X — 1)(a — 2) 
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f k (X) = d k X(A — 1 )(2 — 2)4” c k MX — 1 ) 4~ b k X 4- a k , 

where X u X 2 , 7 3 and X 4 are the four roots of the indicial 
eq. (A12). Using eq. (A10), we have 

7^=4, 7-2 — 3 , 7 3 = 1 , 7 4 = 0. ( A 1 4 ) 

Substituting X x into eq. (A13), w r e obtain the coefficients of the 
first independent solution with an arbitrary constant w 0 to be 
determined. Since all the roots of eq. (A 14) are integers, substi- 
tutions of the rest of the roots into (A13) will not lead to 
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independent solutions. In this ease, the Frobenius technique 
(see Babister 1967) must be invoked for the remaining three 
independent solutions. After the Frobenius operation, we 
complete the four independent solutions for (Al): 


w 0 , Wq, w" and Wq are integral constants to be determined. 
With the vertical component U n known, the horizontal compo- 
nent V n can be derived through the continuity eq. (2), which 
gives 
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The rest of the derivation is straightforward and is thus 
omitted here. 

The correct procedure to deal with the singularities using the 
propagation technique in search of modes is to terminate 
the propagation of eq. ( 5) from the centre somewhere below 
the singularity point r 0 , and then re-initiate two sets ot 
propagation starting at r 0 . One of the propagations goes all 
the way up to the surface, and the other goes down to where 
the original propagation was terminated. The first few terms 
of the solution (A 15) and (A 16) can provide the initial values 
for the latter two propagations. An additional boundary is 
formed at the terminated radius. Flowever, the downward 
propagation just brings an additional set of integral constants 
for that boundary. Hence, the solution of (5) at the surface 
can be uniquely determined. For problems with more than one 
singularity radius, such as the viscosity models C and D in 
Fig. 7, we can proceed following the same principle, that is to 
derive the expansion (A15) and { A 1 6 ) at each singularity point, 
and create a cut-off boundary between each adjacent pair of 
singularity points. 
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